#import gff for anno-TSS
features <- rtracklayer::import("Input/20210217_syne_onlyUnique_withFeat.gff3")
TUs <- rtracklayer::import("Input/Kopf_4091_TUs_combined.gff3")
Coverage data was created using Filtered .bam files with multireads, which were processed by bedtools genome coverage. Files were re-named in the respective history and downloaded as one .zip. In a next step, separate files were combined into one large table using Python code (joinFiles_PSS.py, joinFiles_TSS.py, joinFiles_transcript.py).
Normalization factors are size factors determined DESeq2 using separate .Rmd files (PSS-TSS Analysis.Rmd, TranscriptAnalysis.Rmd).
transcript_factors <- read.table("Input/transcript_sizeFactors.csv", sep=",", header=TRUE, row.names=1)
TSS_factors <- read.table("Input/TSS_sizeFactors.csv", sep=",", header=TRUE, row.names=1)
PSS_factors <- read.table("Input/PSS_sizeFactors.csv", sep=",", header=TRUE, row.names=1)
row.names(TSS_factors) <- c("dWT1", "dWT2", "dWT3", "WT1", "WT2", "WT3", "TV1", "TV2")
row.names(PSS_factors) <- c("dWT1", "dWT2", "dWT3", "WT1", "WT2", "WT3", "TV1", "TV2")
coldata_general <- read.table("Input/colData.csv", sep=",", header=TRUE, row.names=1)
PSS_coverage <- read.table("Input/multireads_PSS_5ends_combined_5sensing.txt", sep="\t", header=TRUE, row.names=1)
PSS_fact <- PSS_factors$factor
if(FALSE){
ddsMat_PSS <- DESeqDataSetFromMatrix(countData = PSS_coverage,
colData = coldata_general,
design = ~ strain)
names(PSS_fact) <- row.names(PSS_factors)
ddsMat_PSS$sizeFactor <- PSS_fact
PSS_norm_counts <- counts(ddsMat_PSS, normalized=TRUE)
grp_File <- data.frame("WT"=apply(PSS_norm_counts[,c("WT1", "WT2", "WT3")],1,mean), "dWT"=apply(PSS_norm_counts[,c("dWT1", "dWT2", "dWT3")],1,mean), "TV"=apply(PSS_norm_counts[,c("TV1", "TV2")],1,mean))
grp_File$seqname <- rep(NA, length(grp_File$WT))
seqnames <- c("BA000022.2", "AP004310.1", "AP004311.1", "AP004312.1", "AP006585.1")
for(i in seqnames){
grp_File[which(grepl(i, row.names(grp_File))),]$seqname <- i
}
grp_File$strand <- rep(NA, length(grp_File$WT))
for(i in c("plus", "minus")){
grp_File[which(grepl(i, row.names(grp_File))),]$strand <- i
}
for(i in seqnames){
for(j in c("plus", "minus")){
tmp <- subset(grp_File, grp_File$seqname==i & grp_File$strand==j)[,c("WT", "dWT", "TV")]
s=""
if(j=="plus"){
s="fw"
}else{
s="rev"
}
write.table(tmp, file=paste("Output/grp/PSS/", i, "_PSS-multireads_WT-dWT-TV_", s, ".grp", sep=""), sep="\t", row.names=FALSE, col.names = FALSE)
}
}
rm(PSS_norm_counts)
}
TSS_coverage <- read.table("Input/multireads_TSS_5ends_combined_5sensing.txt", sep="\t", header=TRUE, row.names=1)
ddsMat_TSS <- DESeqDataSetFromMatrix(countData = TSS_coverage,
colData = coldata_general,
design = ~ strain)
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
TSS_fact <- TSS_factors$factor
names(TSS_fact) <- row.names(TSS_factors)
ddsMat_TSS$sizeFactor <- TSS_fact
TSS_norm_counts <- counts(ddsMat_TSS, normalized=TRUE)
rm(ddsMat_TSS)
if(FALSE){
grp_File <- data.frame("WT"=apply(TSS_norm_counts[,c("WT1", "WT2", "WT3")],1,mean), "dWT"=apply(TSS_norm_counts[,c("dWT1", "dWT2", "dWT3")],1,mean), "TV"=apply(TSS_norm_counts[,c("TV1", "TV2")],1,mean))
grp_File$seqname <- rep(NA, length(grp_File$WT))
seqnames <- c("BA000022.2", "AP004310.1", "AP004311.1", "AP004312.1", "AP006585.1")
for(i in seqnames){
grp_File[which(grepl(i, row.names(grp_File))),]$seqname <- i
}
grp_File$strand <- rep(NA, length(grp_File$WT))
for(i in c("plus", "minus")){
grp_File[which(grepl(i, row.names(grp_File))),]$strand <- i
}
for(i in seqnames){
for(j in c("plus", "minus")){
tmp <- subset(grp_File, grp_File$seqname==i & grp_File$strand==j)[,c("WT", "dWT", "TV")]
s=""
if(j=="plus"){
s="fw"
}else{
s="rev"
}
write.table(tmp, file=paste("Output/grp/TSS/", i, "_TSS-multireads_WT-dWT-TV_", s, ".grp", sep=""), sep="\t", row.names=FALSE, col.names = FALSE)
}
}
}
transcript_coverage <- read.table("Input/multireads_transcript_coverage_combined_5sensing.txt", sep="\t", header=TRUE, row.names=1)
ddsMat_transcript <- DESeqDataSetFromMatrix(countData = transcript_coverage,
colData = coldata_general,
design = ~ strain)
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
transcript_fact <- transcript_factors$factor
names(transcript_fact) <- row.names(transcript_factors)
ddsMat_transcript$sizeFactor <- transcript_fact
transcript_norm_counts <- counts(ddsMat_transcript, normalized=TRUE)
rm(ddsMat_transcript)
rm(transcript_coverage)
if(FALSE){
grp_File <- data.frame("WT"=apply(transcript_norm_counts[,c("WT1", "WT2", "WT3")],1,mean), "dWT"=apply(transcript_norm_counts[,c("dWT1", "dWT2", "dWT3")],1,mean), "TV"=apply(transcript_norm_counts[,c("TV1", "TV2")],1,mean))
grp_File$seqname <- rep(NA, length(grp_File$WT))
seqnames <- c("BA000022.2", "AP004310.1", "AP004311.1", "AP004312.1", "AP006585.1")
for(i in seqnames){
grp_File[which(grepl(i, row.names(grp_File))),]$seqname <- i
}
grp_File$strand <- rep(NA, length(grp_File$WT))
for(i in c("plus", "minus")){
grp_File[which(grepl(i, row.names(grp_File))),]$strand <- i
}
for(i in seqnames){
for(j in c("plus", "minus")){
tmp <- subset(grp_File, grp_File$seqname==i & grp_File$strand==j)[,c("WT", "dWT", "TV")]
s=""
if(j=="plus"){
s="fw"
}else{
s="rev"
}
write.table(tmp, file=paste("Output/grp/transcript/", i, "_transcript-multireads_WT-dWT-TV_", s, ".grp", sep=""), sep="\t", row.names=FALSE, col.names = FALSE)
}
}
}
names(PSS_coverage) <- paste(names(PSS_coverage), "_PSS", sep="")
names(PSS_fact) <- paste(names(PSS_fact), "_PSS", sep="")
names(TSS_coverage) <- paste(names(TSS_coverage), "_TSS", sep="")
names(TSS_fact) <- paste(names(TSS_fact), "_TSS", sep="")
# merge PSS and TSS to merged_raw
merged_raw <- merge(PSS_coverage, TSS_coverage, by="row.names")
row.names(merged_raw) <- merged_raw$Row.names
merged_raw$Row.names <- NULL
rm(TSS_coverage)
When taking into account lowly populated sites, DESeq2 dispersion estimates become over-fitted. Hence, edgeR::filterByExpression() is used as a pre-processing step before further DESeq2 analyses.
group=c(rep("dWT_PSS", 3), rep("WT_PSS", 3), rep("TV_PSS", 2), rep("dWT_TSS", 3), rep("WT_TSS", 3), rep("TV_TSS", 2))
y <- DGEList(counts=merged_raw, group=group)
nrow(y)
## [1] 7894038
keep <- filterByExpr(y)
y <- y[keep, ,keep.lib.size=FALSE]
nrow(y)
## [1] 19944
merged_filtered <- merged_raw[row.names(y$counts),]
nrow(merged_filtered)
## [1] 19944
rm(merged_raw)
coldata <- read.csv("Input/colData_PSS-TSS.csv", row.names=1)
# create DESeq2 data object
ddsMat_PSSTSS <- DESeqDataSetFromMatrix(countData = merged_filtered,
colData = coldata,
design = ~ strain + type)
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
ddsMat_PSSTSS$sizeFactor <- c(PSS_fact, TSS_fact)
ddsMat_PSSTSS <- DESeq(ddsMat_PSSTSS)
## using pre-existing size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
res_PSSTSS <- results(ddsMat_PSSTSS, contrast=c("type", "PSS", "TSS"))
write.csv(res_PSSTSS[order(res_PSSTSS$padj),], file="Output/DESeq2_resultsTables/results_PSS-TSS-comparisons.csv")
plotDispEsts(ddsMat_PSSTSS)
pdf(file="Output/DESeq2_Plots/ddsMat_PSSTSS_DispEsts.pdf", width=4.5, height=4.5)
plotDispEsts(ddsMat_PSSTSS, xlab="Mean of Normalized Counts", ylab="Dispersion")
dev.off()
## png
## 2
p <- PCA_plot_PSSTSS(ddsMat_PSSTSS, "PSS-TSS")
p
ggsave("Output/DESeq2_Plots/ddsMat_PSSTSS_PCA.pdf", plot=p, width=15, height=9, units="cm")
p <- heatmap_plot_PSSTSS(ddsMat_PSSTSS, "PSS-TSS")
p
ggsave("Output/DESeq2_Plots/ddsMat_PSSTSS_heatMap.pdf", plot=p, width=15, height=12, units="cm")
rm(ddsMat_PSSTSS) # to free space
pdf(file="Output/DESeq2_Plots/ddsMat_PSSTSS_hist_log2FC.pdf", width=4.5, height=4.5)
hist(res_PSSTSS$log2FoldChange, breaks=20, main="", xlab="Log2FC(PSS/TSS)", ylab="Frequency")
dev.off()
## png
## 2
hist(res_PSSTSS$log2FoldChange, breaks=20, main="DESeq2 Normalization", xlab="Log2FC PSS/TSS", ylab="Frequency")
# non-normalized
log2_all <- log2(apply(merged_filtered[,1:8],1, mean)/apply(merged_filtered[,9:16],1,mean))
hist(log2_all, breaks=20, main="No normalization", xlab="Log2(PSS/TSS)")
rm(merged_filtered)
p <- MAplot_ggplot(res_PSSTSS, foldchange=0.8, y_axis_label = "Log2 fold-change(PSS/TSS)")
p <- p + scale_colour_manual(values=c("UP"="#4e9879ff", "DOWN"="#9b511fff", "NO"="#d3d3d3ff"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
p
ggsave("Output/DESeq2_Plots/ddsMat_PSSTSS_MAplot.pdf",plot=p, width=15, height=12, units="cm")
res_PSSTSS <- subset(res_PSSTSS, !is.na(res_PSSTSS$padj))
# down: TSS, up: PSS
count_up_down(res_PSSTSS, foldchange=0.8, padjusted=0.05)
## [1] "number of features down: 11620"
## [1] "number of features up: 5424"
p <- volcanoPlot_ggplot(res_PSSTSS, foldchange=0.8, padjusted=0.05)
p
ggsave("Output/DESeq2_Plots/ddsMat_PSSTSS_VolcanoPlot.pdf", plot=p, width=15, height=12, units="cm")
pdf(file="Output/DESeq2_Plots/ddsMat_PSSTSS_pValuePlot.pdf", width=4.5, height=4.5)
pvaluePlot(res_PSSTSS, "Comparison PSS to TSS")
dev.off()
## png
## 2
pvaluePlot(res_PSSTSS, "Comparison PSS to TSS")
# extract which positions are TSS, which are PSS
PSS_positions_unfiltered <- row.names(subset(res_PSSTSS, res_PSSTSS$log2FoldChange>0.8 & res_PSSTSS$padj<0.05))
PSS_nonReduced_Ranges <- create_GRanges_object(base::strsplit(PSS_positions_unfiltered, "-"), res_PSSTSS[PSS_positions_unfiltered,]$baseMean)
TSS_positions_unfiltered <- row.names(subset(res_PSSTSS, res_PSSTSS$log2FoldChange<(-0.8) & res_PSSTSS$padj<0.05))
trans_norm <- as.data.frame(transcript_norm_counts)
rm(transcript_norm_counts)
TSS_norm <- as.data.frame(TSS_norm_counts)
rm(TSS_norm_counts)
TSS_trans_ratio_cutoff <- 0.02
TSS_trans_ratio <- apply(TSS_norm[TSS_positions_unfiltered,],1,mean)/apply(trans_norm[TSS_positions_unfiltered,],1,mean)
TSS_above_cutoff <- subset(TSS_positions_unfiltered, TSS_trans_ratio[TSS_positions_unfiltered]>TSS_trans_ratio_cutoff)
TSS_positions_sorted<- base::strsplit(TSS_above_cutoff, "-")
TSS_nonReduced_Ranges <- create_GRanges_object(TSS_positions_sorted, res_PSSTSS[TSS_above_cutoff,]$baseMean)
rm(trans_norm)
rm(TSS_norm)
PSS_positions_Ranges <- GenomicRanges::reduce(PSS_nonReduced_Ranges)
length(PSS_positions_Ranges)
## [1] 4821
TSS_positions_Ranges <- GenomicRanges::reduce(TSS_nonReduced_Ranges)
length(TSS_positions_Ranges)
## [1] 7406
TSS_positions_Ranges$type <- "TSS"
PSS_positions_Ranges$type <- "PSS"
save_gff(c(TSS_positions_Ranges, PSS_positions_Ranges), "Output/gffs/", "TSS_PSS_multireads")
rm(TSS_positions_Ranges)
rm(PSS_positions_Ranges)
PSS_use <- PSS_coverage[PSS_positions_unfiltered,]
rm(PSS_coverage)
coldata <- read.csv("Input/colData_PSS.csv", row.names=1)
ddsMat_PSS <- DESeqDataSetFromMatrix(countData = PSS_use,
colData = coldata,
design = ~ strain)
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
ddsMat_PSS$sizeFactor <- PSS_fact
ddsMat_PSS <- DESeq(ddsMat_PSS)
## using pre-existing size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
rm(PSS_use)
plotDispEsts(ddsMat_PSS, main="PSS comparison", xlab="Mean of Normalized Counts", ylab="Dispersion")
p <- PCA_plot(ddsMat_PSS, "PSS")
p
ggsave("Output/DESeq2_Plots/ddsMat_PSS_PCA.pdf", plot=p, width=9, height=9, units="cm")
p <- heatmap_plot(ddsMat_PSS, "PSS")
p
ggsave("Output/DESeq2_Plots/ddsMat_PSS_heatMap.pdf", plot=p, width=15, height=12, units="cm")
# extract results
PSS_result_dWT_WT <- results(ddsMat_PSS, contrast=c("strain", "dWT", "WT")) # dWT/WT -> higher in dWT: higher log2FC
write.csv(PSS_result_dWT_WT[order(PSS_result_dWT_WT$padj),], file="Output/DESeq2_resultsTables/results_PSS-dWT-WT.csv")
PSS_result_dWT_TV <- results(ddsMat_PSS, contrast=c("strain", "dWT", "TV")) # dWT/TV -> higher in dWT: higher log2FC
write.csv(PSS_result_dWT_TV[order(PSS_result_dWT_TV$padj),], file="Output/DESeq2_resultsTables/results_PSS-dWT-TV.csv")
PSS_result_WT_TV <- results(ddsMat_PSS, contrast=c("strain", "WT", "TV")) # WT/TV -> higher in WT: higher log2FC
write.csv(PSS_result_WT_TV[order(PSS_result_WT_TV$padj),], file="Output/DESeq2_resultsTables/results_PSS-WT-TV.csv")
rm(ddsMat_PSS)
indices_plus <- which(strand(PSS_nonReduced_Ranges)=="+")
indices_minus <- which(strand(PSS_nonReduced_Ranges)=="-")
PSS_names <- c()
PSS_names[indices_plus] <- paste(seqnames(PSS_nonReduced_Ranges[indices_plus]),start(PSS_nonReduced_Ranges[indices_plus]),"plus",sep="-")
PSS_names[indices_minus] <- paste(seqnames(PSS_nonReduced_Ranges[indices_minus]),start(PSS_nonReduced_Ranges[indices_minus]),"minus",sep="-")
rm(indices_plus)
rm(indices_minus)
PSS_nonReduced_Ranges$names <- PSS_names
PSS_nonReduced_Ranges$featuresOverlap <- rep("", length(PSS_nonReduced_Ranges))
PSS_nonReduced_Ranges$TU_overlap <- rep("", length(PSS_nonReduced_Ranges))
for(i in 1:length(PSS_nonReduced_Ranges)){
PSS_nonReduced_Ranges$featuresOverlap[i] <- paste(features[which(countOverlaps(features,PSS_nonReduced_Ranges[i])>0)]$locus_tag,collapse=",")
PSS_nonReduced_Ranges$TU_overlap[i] <- paste(TUs[which(countOverlaps(TUs,PSS_nonReduced_Ranges[i])>0)]$index,collapse=",")
}
# create tables with annotation: which features and TUs are overlapping with PSS?
df_PSS_nonRed_Ranges <- as.data.frame(PSS_nonReduced_Ranges)
row.names(df_PSS_nonRed_Ranges) <- df_PSS_nonRed_Ranges$names
# PSS result dWT WT
PSS_result_dWT_WT_annot <- rownames_to_column(as.data.frame(PSS_result_dWT_WT[order(PSS_result_dWT_WT$padj),]))
PSS_result_dWT_WT_annot$CDS <- df_PSS_nonRed_Ranges[PSS_result_dWT_WT_annot$rowname,]$featuresOverlap
PSS_result_dWT_WT_annot$TUs <- df_PSS_nonRed_Ranges[PSS_result_dWT_WT_annot$rowname,]$TU_overlap
PSS_result_dWT_WT_annot <- PSS_result_dWT_WT_annot[,c(8,9,1,2:7)]
write_tsv(PSS_result_dWT_WT_annot, "Output/DESeq2_resultsTables/results_PSS_dWT_WT_annotated.tsv")
# PSS result dWT TV
PSS_result_dWT_TV_annot <- rownames_to_column(as.data.frame(PSS_result_dWT_TV[order(PSS_result_dWT_TV$padj),]))
PSS_result_dWT_TV_annot$CDS <- df_PSS_nonRed_Ranges[PSS_result_dWT_TV_annot$rowname,]$featuresOverlap
PSS_result_dWT_TV_annot$TUs <- df_PSS_nonRed_Ranges[PSS_result_dWT_TV_annot$rowname,]$TU_overlap
PSS_result_dWT_TV_annot <- PSS_result_dWT_TV_annot[,c(8,9,1,2:7)]
write_tsv(PSS_result_dWT_TV_annot, "Output/DESeq2_resultsTables/results_PSS_dWT-TV_annotated.tsv")
# PSS result WT, TV
PSS_result_WT_TV_annot <- rownames_to_column(as.data.frame(PSS_result_WT_TV[order(PSS_result_WT_TV$padj),]))
PSS_result_WT_TV_annot$CDS <- df_PSS_nonRed_Ranges[PSS_result_WT_TV_annot$rowname,]$featuresOverlap
PSS_result_WT_TV_annot$TUs <- df_PSS_nonRed_Ranges[PSS_result_WT_TV_annot$rowname,]$TU_overlap
PSS_result_WT_TV_annot <- PSS_result_WT_TV_annot[,c(8,9,1,2:7)]
write_tsv(PSS_result_WT_TV_annot, "Output/DESeq2_resultsTables/results_PSS_WT-TV_annotated.tsv")
Run in respective Directory (output/DESeq2_resultsTables/) python changeAnnotation_DESeq2.py results_PSS_dWT_WT_annotated.tsv results_dWT_WT_annotated-2.tsv python changeAnnotation_DESeq2.py results_PSS_dWT-TV_annotated.tsv results_dWT-TV_annotated-2.tsv python changeAnnotation_DESeq2.py results_PSS_-_WT-TV_annotated.tsv results_PSS-WT-TV_annotated-2.tsv
pdf(file="Output/DESeq2_Plots/ddsMat_PSS-dWT-WT_pValuePlot.pdf", width=4.5, height=4.5)
pvaluePlot(PSS_result_dWT_WT, "PSS dWT WT")
dev.off()
## png
## 2
pdf(file="Output/DESeq2_Plots/ddsMat_PSS-dWT-TV_pValuePlot.pdf", width=4.5, height=4.5)
pvaluePlot(PSS_result_dWT_TV, "PSS dWT dTV")
dev.off()
## png
## 2
pvaluePlot(PSS_result_dWT_WT, "PSS dWT WT")
pvaluePlot(PSS_result_dWT_TV, "PSS dWT dTV")
Filter out PSS below filter level of results(). For these, p.adj is set to NA - hence: remove positions at which p.adj==NA
PSS_result_dWT_WT <- subset(PSS_result_dWT_WT, !is.na(PSS_result_dWT_WT$padj))
PSS_result_dWT_TV <- subset(PSS_result_dWT_TV, !is.na(PSS_result_dWT_TV$padj))
PSS_result_WT_TV <- subset(PSS_result_WT_TV, !is.na(PSS_result_WT_TV$padj))
count_up_down(PSS_result_dWT_WT, foldchange=1, padjusted=0.05)
## [1] "number of features down: 4"
## [1] "number of features up: 67"
p <- volcanoPlot_ggplot(as.data.frame(PSS_result_dWT_WT), foldchange=1, padjusted=0.05) + scale_colour_manual(values=c("DOWN"="#000000ff", "UP"="#005a96ff", "NO"="#d3d3d3ff"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
p
ggsave("Output/DESeq2_Plots/ddsMat_PSS-0h_VolcanoPlot.pdf", plot=p, width=15, height=12, units="cm")
p <- MAplot_ggplot(PSS_result_dWT_WT, foldchange=1.0, y_axis_label = "Log2 fold-change(rne(WT)/WT)") + scale_colour_manual(values=c("DOWN"="#000000ff", "UP"="#005a96ff", "NO"="#d3d3d3ff"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
p
ggsave("Output/DESeq2_Plots/ddsMat_PSS-dWT-WT_MAplot.pdf",plot=p, width=15, height=12, units="cm")
count_up_down(PSS_result_dWT_TV, foldchange=1, padjusted=0.05)
## [1] "number of features down: 175"
## [1] "number of features up: 94"
p <- volcanoPlot_ggplot(as.data.frame(PSS_result_dWT_TV), foldchange=1, padjusted=0.05) + scale_colour_manual(values=c("DOWN"="#e69f00ff", "UP"="#005a96ff", "NO"="#d3d3d3ff"))
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
p
ggsave("Output/DESeq2_Plots/ddsMat_PSS-dWT-TV_VolcanoPlot.pdf", plot=p, width=15, height=12, units="cm")
p <- MAplot_ggplot(PSS_result_dWT_TV, foldchange=1.0, y_axis_label = "Log2 fold-change(rne(WT)/rne(5p))") + scale_colour_manual(values=c("DOWN"="#e69f00ff", "UP"="#005a96ff", "NO"="#d3d3d3ff"))# + ylim(-5,+5)
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
p
ggsave("Output/DESeq2_Plots/ddsMat_PSS-1h_MAplot.pdf",plot=p, width=15, height=12, units="cm")
count_up_down(PSS_result_WT_TV, foldchange = 1, padjusted = 0.05)
## [1] "number of features down: 237"
## [1] "number of features up: 111"
p <- MAplot_ggplot(PSS_result_WT_TV, foldchange=1.0, y_axis_label = "Log2 fold-change(WT/rne(5p))") + scale_colour_manual(values=c("DOWN"="#e69f00ff", "UP"="#000000ff", "NO"="#d3d3d3ff"))# + ylim(-5,+5)
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
p
ggsave("Output/DESeq2_Plots/ddsMat_PSS-WT-TV_MAplot.pdf",plot=p, width=15, height=12, units="cm")
## R version 4.0.4 (2021-02-15)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Debian GNU/Linux 11 (bullseye)
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
## LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
##
## locale:
## [1] LC_CTYPE=de_DE.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=de_DE.UTF-8 LC_COLLATE=de_DE.UTF-8
## [5] LC_MONETARY=de_DE.UTF-8 LC_MESSAGES=de_DE.UTF-8
## [7] LC_PAPER=de_DE.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=de_DE.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] viridis_0.5.1 viridisLite_0.3.0
## [3] rtracklayer_1.50.0 Biostrings_2.58.0
## [5] XVector_0.30.0 vsn_3.58.0
## [7] RColorBrewer_1.1-2 pheatmap_1.0.12
## [9] ggpubr_0.4.0 ggrepel_0.9.1
## [11] edgeR_3.32.1 limma_3.46.0
## [13] DESeq2_1.30.1 SummarizedExperiment_1.20.0
## [15] Biobase_2.50.0 MatrixGenerics_1.2.1
## [17] matrixStats_0.58.0 GenomicRanges_1.42.0
## [19] GenomeInfoDb_1.26.2 IRanges_2.24.1
## [21] S4Vectors_0.28.1 BiocGenerics_0.36.0
## [23] forcats_0.5.1 stringr_1.4.0
## [25] dplyr_1.0.5 purrr_0.3.4
## [27] readr_1.4.0 tidyr_1.1.3
## [29] tibble_3.1.0 ggplot2_3.3.3
## [31] tidyverse_1.3.0 knitr_1.31
##
## loaded via a namespace (and not attached):
## [1] colorspace_2.0-0 ggsignif_0.6.1 ellipsis_0.3.1
## [4] rio_0.5.26 fs_1.5.0 rstudioapi_0.13
## [7] farver_2.1.0 affyio_1.60.0 bit64_4.0.5
## [10] AnnotationDbi_1.52.0 fansi_0.4.2 lubridate_1.7.10
## [13] xml2_1.3.2 splines_4.0.4 cachem_1.0.4
## [16] geneplotter_1.68.0 jsonlite_1.7.2 Rsamtools_2.6.0
## [19] broom_0.7.5 annotate_1.68.0 dbplyr_2.1.0
## [22] BiocManager_1.30.10 compiler_4.0.4 httr_1.4.2
## [25] backports_1.2.1 assertthat_0.2.1 Matrix_1.3-2
## [28] fastmap_1.1.0 cli_2.3.1 htmltools_0.5.1.1
## [31] tools_4.0.4 affy_1.68.0 gtable_0.3.0
## [34] glue_1.4.2 GenomeInfoDbData_1.2.4 Rcpp_1.0.6
## [37] carData_3.0-4 cellranger_1.1.0 jquerylib_0.1.3
## [40] vctrs_0.3.6 preprocessCore_1.52.1 xfun_0.22
## [43] openxlsx_4.2.3 rvest_1.0.0 lifecycle_1.0.0
## [46] rstatix_0.7.0 XML_3.99-0.5 zlibbioc_1.36.0
## [49] scales_1.1.1 hms_1.0.0 yaml_2.2.1
## [52] curl_4.3 gridExtra_2.3 memoise_2.0.0
## [55] sass_0.3.1 stringi_1.5.3 RSQLite_2.2.3
## [58] highr_0.8 genefilter_1.72.1 zip_2.1.1
## [61] BiocParallel_1.24.1 rlang_0.4.10 pkgconfig_2.0.3
## [64] bitops_1.0-6 evaluate_0.14 lattice_0.20-41
## [67] labeling_0.4.2 GenomicAlignments_1.26.0 bit_4.0.4
## [70] tidyselect_1.1.0 magrittr_2.0.1 R6_2.5.0
## [73] generics_0.1.0 DelayedArray_0.16.2 DBI_1.1.1
## [76] pillar_1.5.1 haven_2.3.1 foreign_0.8-81
## [79] withr_2.4.1 abind_1.4-5 survival_3.2-7
## [82] RCurl_1.98-1.2 modelr_0.1.8 crayon_1.4.1
## [85] car_3.0-10 utf8_1.1.4 rmarkdown_2.7
## [88] locfit_1.5-9.4 grid_4.0.4 readxl_1.3.1
## [91] data.table_1.14.0 blob_1.2.1 reprex_1.0.0
## [94] digest_0.6.27 xtable_1.8-4 munsell_0.5.0
## [97] bslib_0.2.4